#### Start of replication of cross-national analyses ####

# Load packages
library(readr)
library(lfe)
library(nlme)
library(dplyr)
library(stargazer)
library(car)
library(plm)
library(marginaleffects)
library(giscoR)
library(countrycode)
library(ggplot2)
library(sf)
library(ivDiag)


# Clear environment
rm(list = ls())

# @USER: PLEASE CHANGE TO YOUR WORKING DIRECTORY
#setwd("C:/Users/HSmidt/Dropbox/Projekte/Paper_05 - CSO repression and Human Rights complaints/03_Code_and_Data/")
setwd("C:/Users/chrstein/Dropbox/CSO repression and Human Rights complaints/03_Code_and_Data/")
 
# Read data for analysis
load("10_replication_files/data_for_analysis_crossnat.RData")


#### Check media freedom and democracy in India ####

# Number communications with and without signature of rapporteur on human rights defenders
sum(data$mandate_hrdefenders[data$cowc=="IND"]) / ( sum(data$mandate_nohrdefenders[data$cowc=="IND"]) + sum(data$mandate_hrdefenders[data$cowc=="IND"]) )

# Number communications with and without signature of rapporteur on human rights defenders
sum(data$mandate_hrdefenders) / ( sum(data$mandate_nohrdefenders) + sum(data$mandate_hrdefenders) )

# Media freedom 
mediafree <- data %>% dplyr::select(year, country, v2mecenefm) %>% 
  filter(year > 2009) %>% filter(year < 2023)
mean(mediafree$v2mecenefm)
mediafree_grouped <- mediafree %>% group_by(country) %>% 
  summarize(avg = mean(v2mecenefm, na.rm = T)) %>% arrange(avg)
tercile <- quantile(mediafree_grouped$avg, probs = seq(0, 1, by = 1/3))
mediafree_grouped$tercile <- cut(mediafree_grouped$avg, breaks = tercile, labels = F)
quintile <- quantile(mediafree_grouped$avg, probs = seq(0, 1, by = 1/5))
mediafree_grouped$quintile <- cut(mediafree_grouped$avg, breaks = quintile, labels = F)

# Democracy 
democracy <- data %>% dplyr::select(year, country, v2x_polyarchy) %>% 
  filter(year > 2009) %>% filter(year < 2023)
mean(democracy$v2x_polyarchy)
democracy_grouped <- democracy %>% group_by(country) %>% 
  summarize(avg = mean(v2x_polyarchy, na.rm = T)) %>% arrange(avg)
tercile <- quantile(democracy_grouped$avg, probs = seq(0, 1, by = 1/3))
democracy_grouped$tercile <- cut(mediafree_grouped$avg, breaks = tercile, labels = F)
quintile <- quantile(democracy_grouped$avg, probs = seq(0, 1, by = 1/5))
democracy_grouped$quintile <- cut(democracy_grouped$avg, breaks = quintile, labels = F)


#*##########################################################*#  
#### Figure 1 (main ms.) and Figures A1 and A2 (App. C) #####
#*##########################################################*#

##### Map Figure 1, main ms. early average of number of UNSP communications from 2010 to 2021 #####

# Get coordinates for maps
countries <- gisco_get_countries(year = "2016") # does only work for this year
countries <- subset(countries, CNTR_NAME != "Antarctica")
countries_small <- dplyr::select(countries, ISO3_CODE, geometry)

# Count the number of UNSP Communications per country
country_based_comm <- data %>% 
  filter(year > 2009) %>%
  filter(year < 2022) %>% 
  group_by(country) %>% 
  summarise(commcount_individual = mean(commcount_individual, na.rm=T) ) %>% 
  arrange(desc(commcount_individual)) %>%
  ungroup()

# Get coordinates and merge
country_based_comm$ISO3_CODE <- countrycode(country_based_comm$country, origin = "country.name", destination = "iso3c")
country_based_comm <- merge(country_based_comm, countries_small, by = "ISO3_CODE", all.x = T, all.y = T) 

# Recode NAs
country_based_comm$commcount_individual <- ifelse(is.na(country_based_comm$commcount_individual), 0, country_based_comm$commcount_individual)

# Create World Map of UNSP Communications (Figure 5)
ragg::agg_png("10_replication_files/Outputs/Figures/Fig1_mainms_map_unspcomm.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(data = country_based_comm$geometry) + 
  geom_sf(aes(fill = country_based_comm$commcount_individual), size = 0.02) +
  scale_fill_viridis_c(option = "A", alpha = 1, direction = -1, name = "Average number of \nUNSP communications") +
  theme(plot.title = element_text(hjust = 0.5)) +  
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
dev.off()
 

##### Map A1, App C Country averages of CSO repression (V-Dem), 2011–2022 #####
 
# Make average repression
cso_repression_avg <- data %>% 
  group_by(country) %>% 
  summarise(cso_repress_avg = mean(v2csreprss,na.rm =T)) %>% 
  filter(cso_repress_avg != "NaN") %>% 
  arrange(desc(cso_repress_avg)) %>%
  ungroup()

# Get Coordinates and Merge
cso_repression_avg$ISO3_CODE <- countrycode(cso_repression_avg$country, origin = "country.name", destination = "iso3c")
cso_repression_avg <- merge(cso_repression_avg, countries_small, by = "ISO3_CODE", all.x = T, all.y = T) 

# Create World Map of CSO Repression
ragg::agg_png("10_replication_files/Outputs/Figures/FigA1_AppC_map_csorepression.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(data = cso_repression_avg$geometry) + 
  geom_sf(aes(fill = cso_repression_avg$cso_repress_avg), size = 0.02) +
  scale_fill_viridis_c(option = "A", alpha = 1, direction = -1, name = "Average value \nCSO repression") +
  theme(plot.title = element_text(hjust = 0.5)) +  
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
dev.off()

##### Map A2, App C Country averages of reprisal for UN cooperation (original data), 2011–2022 #####

# Count reported reprisals for UN cooperation
country_based_reprisals <- data %>% 
  group_by(country) %>% 
  summarise(unreprisals = sum(unreprisal_acc)) %>% 
  arrange(desc(unreprisals)) %>%
  ungroup()

# Merge with coordinates
country_based_reprisals$ISO3_CODE <- countrycode(country_based_reprisals$country, origin = "country.name", destination = "iso3c")
country_based_reprisals <- merge(country_based_reprisals, countries_small, by = "ISO3_CODE", all.x = T, all.y = T) 

# Recode NAs
country_based_reprisals$unreprisals <- ifelse(is.na(country_based_reprisals$unreprisals), 0, country_based_reprisals$unreprisals)

# Create World Map of Reprisals 
ragg::agg_png("10_replication_files/Outputs/Figures/FigA2_AppC_map_unreprisals.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(data = country_based_reprisals$geometry) + 
  geom_sf(aes(fill = country_based_reprisals$unreprisals), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of years with \n reprisals") +
  theme(plot.title = element_text(hjust = 0.5)) +  
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
dev.off()


#*###############################*#
#### Define control variables #####
#*###############################*#


# Test for multicollinearity #
mtest <- lm(v2csreprss ~  commcount_individual_lag1 + v2peedueq_lag1 + conflict_lag + v2mecenefm_lag1 + v2juhcind_lag1 + log_gdpcap_lag1, data=data)
summary(mtest)
vif(mtest)

# Create vector with controls for cross-country models  
controls <- paste0("+ v2peedueq_lag1 + conflict_lag + v2mecenefm_lag1 + v2juhcind_lag1 + log_gdpcap_lag1")


#*###########################################*#
#### First Stage Model in Main Body of Ms. ####
#*###########################################*#

# Main model
f4 <- as.formula( paste0("v2csreprss ~  1 ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))
m4 <- felm(formula = f4, data = data)

# First stage model
first <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
summary(first)

fstat1 <- round(summary(first)$fstat, 3)
fstat1

custom_significance <- c("\\dag", "*", "**", "***")
star_first_full <- stargazer(first
                             , title = "First stage model"
                             , digits=3
                             , font.size="scriptsize"
                             , omit.stat = c("n", "ser")
                             , star.char = custom_significance
                             , star.cutoffs = c(0.1, 0.05, 0.01, 0.001) 
                             , add.lines = list(c("Number of observations"
                                                , length(first$fitted.values)) ,
                                                c("Country fixed effects", "Yes"),
                                                c("Time-varying controls", "Yes"),
                                                c("Lagged depend. var.", "No"),
                                                c("First stage F-statistic", fstat1))
                             , omit = c("v2peedueq_lag1",
                                        "conflict_lag",
                                        "v2mecenefm_lag1",
                                        "v2juhcind_lag1",
                                        "log_gdpcap_lag1",
                                        "spvisits_lag1")
                             , covariate.labels = c("New UNSP experts (lag 1 yr)"
                                                    , "Number of UNSP experts  (lag 1 yr)"
                                                    )
                             , model.names=F, model.numbers=F
                             , dep.var.labels = "DV: UNSP communications"
                             , label = "tab:first_stage_country"
                             , table.placement="H" 
                             , column.sep.width = "1pt"
                             , align=F )

note.latex <- "\\multicolumn{2}{l} {\\parbox[t]{10cm}{ \\textit{Notes:} The models are estimated with ordinary least squares regression. The control variables are education (lag 1 yr), conflict occurrence (lag 1 yr), freedom of media (lag 1 yr), judicial independence (lag 1 yr), logged GDP per capita (lag 1 yr), and UNSP country visits (lag 1 yr). Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_first_full[grepl("Note",star_first_full)] <- note.latex
cat(star_first_full, sep = "\n")
write(star_first_full, "./10_replication_files/Outputs/Tables/Table1_mainms_first_stage_country.tex")




#*################################*#
#### Models in Main Body of Ms. ####
#*################################*#

#### Main models of repression of CSOs ####

f0 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  1 ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))
 
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# First stage (save for analysis_V3_first_stage.R to make the table with all first stage regression models)
first <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )

# Summary statistics
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  


##### Table A.1 for Appendix C #####

custom_significance <- c("\\dag", "*", "**", "***")
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
write(star_full, "./10_replication_files/Outputs/Tables/TableA1_AppC_v2csreprss_full.tex")


##### Figure 2 (main ms) #####

# Marginal effect of communications
effect_sizes_csos <- matrix(NA, nrow=5, ncol=5)
coef_m0 <- MASS::mvrnorm(1000, coefficients(m0), vcov(m0))[,2] / sd( model.frame(m0)[,c("v2csreprss")] )
effect_sizes_csos[1,] <- quantile(coef_m0, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m1 <- MASS::mvrnorm(1000, coefficients(m1), vcov(m1))[,1] / sd( model.frame(m1)[,c("v2csreprss")] )
effect_sizes_csos[2,] <- quantile(coef_m1, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m2 <- MASS::mvrnorm(1000, coefficients(m2), vcov(m2))[,1] / sd( model.frame(m2)[,c("v2csreprss")] )
effect_sizes_csos[3,] <- quantile(coef_m2, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m3 <- MASS::mvrnorm(1000, coefficients(m3), vcov(m3))[,1] / sd( model.frame(m2)[,c("v2csreprss")] )
effect_sizes_csos[4,] <- quantile(coef_m3, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m4 <- MASS::mvrnorm(1000, coefficients(m4), vcov(m4))[,7] / sd( model.frame(m4)[,c("v2csreprss")] )
effect_sizes_csos[5,] <- quantile(coef_m4, c(0.025, 0.05, 0.5, 0.95, 0.975))

# Compare effect of communication to effect of armed conflict
coefficients(m4)[7] / sd( model.frame(m4)[,c("v2csreprss")] ) 
coefficients(m4)[2] / sd( model.frame(m4)[,c("v2csreprss")] ) 

# Graph 
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig2_mainms_marg_effect_crossnat_csorepress.png"
               , width = 297, height = 210, units = "mm", res = 600)
par(mar = c(6.1, 0.1, 2.1, 4.1),mgp = c(3.5, 1, 0), xpd=T)
plot( y = c(1,2,3,4,5)
      , x = rev(effect_sizes_csos[,3])*100
      , pch= rev(c(16,15,17,18,20))
      , cex=1.5
      , type="p"
      , axes= F
      , xlab = "Effect of one additional UNSP communication (in % standard deviation of CSO repression)"
      , cex.lab = 1.5
      , ylab = ""
      , ylim = c(0.5, 5.5), xlim = c(-15,15)  )
arrows( x0 = rev(effect_sizes_csos[,1]*100), x1 = rev(effect_sizes_csos[,5]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=1)
arrows( x0 = rev(effect_sizes_csos[,2]*100), x1 = rev(effect_sizes_csos[,4]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=2)
lines(y=c(0.5,5.5), x = c(0,0), lty="dashed")
axis(1, at = seq(-10, 15, by=5), labels = seq(-10, 15, by=5), cex.axis=1.2)
text(x = rep(-10,5), y=c(1,2,3,4,5)
     , labels=c("With CSO FEs, \ncontrols, and \ninstrument"
                , "With CSO FEs, \ncontrols, and \ninstrumented lagged DV"
                , "With CSO FEs and \ncontrols"
                , "With CSO FEs"
                , "Bivariate regression" ), cex = 1.5, pos = 4)
text(x= -10, y=5.5, labels = "Models:", font=2, cex=1.5, pos=4)
dev.off()


##### Sensitivity check for IV regression #####

Y <- "v2csreprss" # Y: outcome of interest
D <-"commcount_individual_lag1" # D: endogenous treatment
Z <- c("new_entries_sum_lag1", "Rapporteurs_active_number_lag1") # Z: instrumental variable
control_vars <- c("v2peedueq_lag1", "conflict_lag", "v2mecenefm_lag1", "v2juhcind_lag1", "log_gdpcap_lag1","spvisits_lag1") 
FE <- c("cowc") # FE
cl <- c("cowc") # clusters
g <- ivDiag(data = data, Y=Y, D = D, Z = Z, controls = control_vars, FE= FE, cl = cl)
g$est_ols
g$est_2sls
g$AR
g$F_stat
g$rho
g$est_fs
g$est_rf



#*#######################################################*#
#### Main models of reprisals against UN collaborators ####
#*#######################################################*#

##### Table A.2 for Appendix C #####

# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1 | cowc | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
summary(first)
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  



# Table for Appendix
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications (lag 1 yr,  instr.)")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
write(star_full, "./10_replication_files/Outputs/Tables/TableA2_AppC_unreprisal_acc_full.tex")



##### Figure 3 (main ms) #####

# Effects of communications
effect_sizes_csos <- matrix(NA, nrow=5, ncol=5)
coef_m0 <- MASS::mvrnorm(1000, coefficients(m0), vcov(m0))[,1] / sd( model.frame(m0)[,c("unreprisal_acc")] )
effect_sizes_csos[1,] <- quantile(coef_m0, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m1 <- MASS::mvrnorm(1000, coefficients(m1), vcov(m1))[,1] / sd( model.frame(m1)[,c("unreprisal_acc")] )
effect_sizes_csos[2,] <- quantile(coef_m1, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m2 <- MASS::mvrnorm(1000, coefficients(m2), vcov(m2))[,1] / sd( model.frame(m2)[,c("unreprisal_acc")] )
effect_sizes_csos[3,] <- quantile(coef_m2, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m3 <- MASS::mvrnorm(1000, coefficients(m3), vcov(m3))[,1] / sd( model.frame(m2)[,c("unreprisal_acc")] )
effect_sizes_csos[4,] <- quantile(coef_m3, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m4 <- MASS::mvrnorm(1000, coefficients(m4), vcov(m4))[,7] / sd( model.frame(m4)[,c("unreprisal_acc")] )
effect_sizes_csos[5,] <- quantile(coef_m4, c(0.025, 0.05, 0.5, 0.95, 0.975))

# Compare effect of communication to effect of armed conflict
( coefficients(m4)[7] / sd( model.frame(m4)[,c("unreprisal_acc")] ) ) 
( coefficients(m4)[2] / sd( model.frame(m4)[,c("unreprisal_acc")] ) )

# Graph
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig3_mainms_marg_effect_crossnat_reprisal.png"
               , width = 297, height = 210, units = "mm", res = 600)

par(mar = c(5.1, 4.1, 2.1, 5.1), mgp = c(3.5, 1, 0), xpd=T)
plot( y = c(1,2,3,4,5)
      , x = rev(effect_sizes_csos[,3]*100)
      , pch= rev(c(16,15,17,18,20))
      , cex=1.5
      , type="p"
      , axes= F
      , cex.lab = 1.5
      , xlab = "Effect of one additional UNSP communication (in % standard deviation of CSO repression)"
      , ylab = ""
      , ylim = c(0.5, 5.5), xlim = c(-20,30) )
arrows( x0 = rev(effect_sizes_csos[,1]*100), x1 = rev(effect_sizes_csos[,5]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=1)
arrows( x0 = rev(effect_sizes_csos[,2]*100), x1 = rev(effect_sizes_csos[,4]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=2)
lines(y=c(0.5,5.5), x = c(0,0), lty="dashed")
axis(1, at = seq(-20, 30, by=10), labels = seq(-20, 30, by=10), cex.axis=1.2)
text(x = rep(-20,5), y=c(1,2,3,4,5)
     , labels=c("With country FEs, \ncontrols, and \ninstrument"
                , "With country FEs, \ncontrols, and \ninstrumented lagged DV"
                , "With country FEs and \ncontrols"
                , "With country FEs"
                , "Bivariate regression" ), cex = 1.5, pos = 4)
text(x= -20, y=5.5, labels = "Models:", font=2, cex=1.5, pos=4)

dev.off()


##### Sensitivity check for IV #####
Y <- "unreprisal_acc" # Y: outcome of interest
g <- ivDiag(data = data, Y=Y, D = D, Z = Z, controls = control_vars, FE= FE, cl = cl)
g$est_ols
g$est_2sls
g$AR
g$F_stat
g$rho
g$est_fs
g$est_rf




#*######################################*#
          #### APPENDIX ####
#*######################################*#


#### Binary UNSP communications ####

f0 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1_bi | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1_bi | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1_bi", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1_bi", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1_bi ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# Summary statistics
first <- felm(as.formula( paste0("commcount_individual_lag1_bi ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.5 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and occurrence of at least one UNSP communication"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications binary (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications binary (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full_binary"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_binary.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA5_AppC_v2csreprss_full_binary.tex")



# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1_bi  | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1_bi | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1_bi", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1_bi", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (commcount_individual_lag1_bi ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("commcount_individual_lag1_bi ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.6 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and occurrence of at least one UNSP communication"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications binary (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications binary (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full_binary"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_binary.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA6_AppC_unreprisal_acc_full_binary.tex")



#### Only UNSP communications regarding physical integrity violations ####

f0 <- as.formula( paste0("v2csreprss ~  mandate_physicalintegr_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  mandate_physicalintegr_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  mandate_physicalintegr_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  mandate_physicalintegr_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  ", controls , " + spvisits_lag1 | cowc | (mandate_physicalintegr_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# Summary statistics
first <- felm(as.formula( paste0("mandate_physicalintegr_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.7 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP communications regarding physical integrity rights"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications on phys. integr (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications on phys. integr (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full_physint"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_physint.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA7_AppC_v2csreprss_full_physint.tex")



# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ mandate_physicalintegr_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ mandate_physicalintegr_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ mandate_physicalintegr_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ mandate_physicalintegr_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (mandate_physicalintegr_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("mandate_physicalintegr_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.8 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP communications regarding physical integrity rights"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications on phys. integr (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications on phys. integr (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full_physint"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_physint.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA8_AppC_unreprisal_acc_full_physint.tex")



#### Only UNSP communications that are NOT signed by rapporteurs on human rights defenders  ####
sum(data$mandate_hrdefenders) / sum(data$commcount_individual)
sum(data$mandate_nohrdefenders) / sum(data$commcount_individual)

f0 <- as.formula( paste0("v2csreprss ~  mandate_nohrdefenders_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  mandate_nohrdefenders_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  mandate_nohrdefenders_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  mandate_nohrdefenders_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  ", controls , " + spvisits_lag1 | cowc | (mandate_nohrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# Summary statistics
first <- felm(as.formula( paste0("mandate_nohrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.11 for Appendix C #####
custom_significance <- c("\\dag", "*", "**", "***")
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP communications not regarding human rights defenders"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications not on hr defenders (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications not on hr defenders (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full_nohrdefenders"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_nohrdefenders.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA11_AppC_v2csreprss_full_nohrdefenders.tex")



# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ mandate_nohrdefenders_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ mandate_nohrdefenders_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ mandate_nohrdefenders_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ mandate_nohrdefenders_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (mandate_nohrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("mandate_nohrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.12 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP communications not regarding human rights defenders"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications not on hr defenders (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications not on hr defenders (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full_nohrdefenders"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_nohrdefenders.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA12_AppC_unreprisal_acc_full_nohrdefenders.tex")



#### Only UNSP communications that ARE signed by rapporteurs on human rights defenders  ####
sum(data$mandate_hrdefenders) / sum(data$commcount_individual)
sum(data$mandate_nohrdefenders) / sum(data$commcount_individual)

f0 <- as.formula( paste0("v2csreprss ~  mandate_hrdefenders_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  mandate_hrdefenders_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  mandate_hrdefenders_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  mandate_hrdefenders_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  ", controls , " + spvisits_lag1 | cowc | (mandate_hrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# Summary statistics
first <- felm(as.formula( paste0("mandate_hrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

# Table not displayed
custom_significance <- c("\\dag", "*", "**", "***")
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP communications regarding human rights defenders"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications on hr defenders (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged gdp per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications on hr defenders (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full_hrdefenders"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_hrdefenders.tex")


# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ mandate_hrdefenders_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ mandate_hrdefenders_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ mandate_hrdefenders_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ mandate_hrdefenders_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (mandate_hrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("mandate_hrdefenders_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

# Table not displayed
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP communications regarding human rights defenders"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications on hr defenders (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged gdp per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications on hr defenders (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full_hrdefenders"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_hrdefenders.tex")



#### Only UNSP communications that are retrospective, allegation letters ####

f0 <- as.formula( paste0("v2csreprss ~  allegation_letters_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  allegation_letters_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  allegation_letters_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  allegation_letters_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  ", controls , " + spvisits_lag1 | cowc | (allegation_letters_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# Summary statistics
first <- felm(as.formula( paste0("allegation_letters_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.9 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP allegation letters"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("Allegation letters (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged gdp per capita (lag 1 yr)"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "Allegation letters (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full_allegation"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_allegation.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA9_AppC_v2csreprss_full_allegation.tex")


# Formulas
f0 <- as.formula( paste0("unreprisal_acc ~ allegation_letters_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ allegation_letters_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ allegation_letters_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ allegation_letters_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (allegation_letters_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("allegation_letters_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.10 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP allegation letters"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("Allegation letters (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged gdp per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "Allegation letters (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full_allegation"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_allegation.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA10_AppC_unreprisal_acc_full_allegation.tex")



#### Subset Analyses by Polyarchy Median Split ####

## Repression of CSOs 

f3 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  1 ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m3d <- pgmm(formula = f3, data = subset(data, v2x_polyarchy>=median(v2x_polyarchy)), index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4d <- felm(formula = f4, data = subset(data, v2x_polyarchy>=median(v2x_polyarchy)) )

m3a <- pgmm(formula = f3, data = subset(data, v2x_polyarchy<median(v2x_polyarchy)), index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4a <- felm(formula = f4, data = subset(data, v2x_polyarchy<median(v2x_polyarchy)) )

# Summary statistics
first_d <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4d) )
fstat_d <- round(summary(first_d)$fstat, 3)
fstat_d
first_a <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4a) )
fstat_a <- round(summary(first_a)$fstat, 3)
fstat_a

##### Table A.13 for Appendix C #####
star_reduced <- stargazer(m3d, m4d, m3a, m4a
                          , type="latex"
                          , title = "Repression of CSOs and UNSP communications in sub-samples"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m4d$fitted.values)
                                               , length(m4d$fitted.values)
                                               , length(m4a$fitted.values)
                                               , length(m4a$fitted.values) ) ,
                                             c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "Yes", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "Yes", "No","Yes", "No"),
                                             c("First stage F-statistic", "NA", fstat_d, "NA",  fstat_a))
                          , keep =  "commcount_individual_lag1"
                          , covariate.labels = c("UNSP communications (lag 1 yr)"
                                                 , "UNSP communications (lag 1 yr,  instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)" )
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:v2csreprss_reduced_sub"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The model in columns 2 and 4 are estimated with ordinary least squares regression, while the model in columns 1 and 3 are estimated using generalized methods of moments. The lagged dependent variable in the models in column 1 and 3 is instrumented with deeper lags. In all models, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year. In the models in columns 2 and 4, we additional control for UNSP country visits. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./outputs/v2csreprss_reduced_sub.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA13_AppC_v2csreprss_reduced_sub.tex")


# Table not displayed
star_full <- stargazer(m3d, m4d, m3a, m4a
                          , type="latex"
                          , title = "Repression of CSOs and UNSP communications in sub-samples"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m4d$fitted.values)
                                               , length(m4d$fitted.values)
                                               , length(m4a$fitted.values)
                                               , length(m4a$fitted.values) ) ,
                                             c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "Yes", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "Yes", "No","Yes", "No"),
                                             c("First stage F-statistic", "NA", fstat_d, "NA",  fstat_a))
                          , covariate.labels = c("UNSP communications (lag 1 yr)"
                                                 , "Education (lag 1 yr)"
                                                 , "Conflict (lag 1 yr)"
                                                 , "Freedom of media (lag 1 yr)"
                                                 , "Judicial independence (lag 1 yr)"
                                                 , "Logged GDP per capita (lag 1 yr)"
                                                 , "Repression of CSOs (lag 1 yr, instr.)"
                                                 , "UNSP country visits (lag 1 yr)"
                                                 , "UNSP communications (lag 1 yr,  instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)" )
                          , dep.var.labels.include = F
                          , dep.var.caption = "Dependent variable: Repression of CSOs (V-Dem)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:v2csreprss_full_sub"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{4}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The model in columns 2 and 4 are estimated with ordinary least squares regression, while the model in columns 1 and 3 are estimated using generalized methods of moments. The lagged dependent variable in the models in column 1 and 3 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_sub.tex")


## Reprisals for cooperation with UN

f3 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~  1 ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m3d <- pgmm(formula = f3, data = subset(data, v2x_polyarchy>=median(v2x_polyarchy)), index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4d <- felm(formula = f4, data = subset(data, v2x_polyarchy>=median(v2x_polyarchy)) )

m3a <- pgmm(formula = f3, data = subset(data, v2x_polyarchy<median(v2x_polyarchy)), index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4a <- felm(formula = f4, data = subset(data, v2x_polyarchy<median(v2x_polyarchy)) )

# Summary statistics
first_d <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4d) )
fstat_d <- round(summary(first_d)$fstat, 3)
fstat_d
first_a <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4a) )
fstat_a <- round(summary(first_a)$fstat, 3)
fstat_a

##### Table A.14 for Appendix C #####
star_reduced <- stargazer(m3d, m4d, m3a, m4a
                          , type="latex"
                          , title = "Reprisal for cooperation with UN and UNSP communications in sub-samples"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m4d$fitted.values)
                                               , length(m4d$fitted.values)
                                               , length(m4a$fitted.values)
                                               , length(m4a$fitted.values) ) ,
                                             c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "Yes", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "Yes", "No","Yes", "No"),
                                             c("First stage F-statistic", "NA", fstat_d, "NA",  fstat_a))
                          , keep =  "commcount_individual_lag1"
                          , covariate.labels = c("Communications (lag 1 yr)"
                                                 , "Communications (lag 1 yr,  instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)" )
                          , dep.var.labels.include = F
                          , dep.var.caption = "Dependent variable: Reprisal (original data)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:unreprisal_acc_reduce_sub"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{15.8cm}{ \\textit{Notes:} The model in columns 2 and 4 are estimated with ordinary least squares regression, while the model in columns 1 and 3 are estimated using generalized methods of moments. The lagged dependent variable in the models in column 1 and 3 is instrumented with deeper lags. In all models, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year. In the models in columns 2 and 4, we additional control for UNSP country visits. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./outputs/unreprisal_acc_reduce_sub.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA14_AppC_unreprisal_acc_reduce_sub.tex")


# Table not displayed
star_full <- stargazer(m3d, m4d, m3a, m4a
                          , type="latex"
                          , title = "Reprisal for cooperation with UN and UNSP communications in sub-samples"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m4d$fitted.values)
                                               , length(m4d$fitted.values)
                                               , length(m4a$fitted.values)
                                               , length(m4a$fitted.values) ) ,
                                             c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "Yes", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "Yes", "No","Yes", "No"),
                                             c("First stage F-statistic", "NA", fstat_d, "NA",  fstat_a))
                          , covariate.labels = c("UNSP communications (lag 1 yr)"
                                                 , "Education (lag 1 yr)"
                                                 , "Conflict (lag 1 yr)"
                                                 , "Freedom of media (lag 1 yr)"
                                                 , "Judicial independence (lag 1 yr)"
                                                 , "Logged GDP per capita (lag 1 yr)"
                                                 , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                                 , "UNSP country visits (lag 1 yr)"
                                                 , "UNSP communications (lag 1 yr,  instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)" )
                          , dep.var.labels.include = F
                          , dep.var.caption = "Dependent variable: Reprisal (original data)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:unreprisal_acc_full_sub"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{4}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The model in columns 2 and 4 are estimated with ordinary least squares regression, while the model in columns 1 and 3 are estimated using generalized methods of moments. The lagged dependent variable in the models in column 1 and 3 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_sub.tex")





#### Interaction of Polyarchy and Communications ####

# CSO repression
f0 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 * v2x_polyarchy_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 * v2x_polyarchy_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 * v2x_polyarchy_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 * v2x_polyarchy_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')


##### Table A.15 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3
                       , type="latex"
                       , title = "Repression of CSOs and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values) ) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes"))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Polyarchy (1 yr lag)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "CSO repression (lag 1 yr, instr.)"
                                              , "UNSP communications * Polyarchy"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:v2csreprss_full_interaction"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{13cm}{ \\textit{Notes:} The models in columns 1, 2 and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_reduced_interaction.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA15_AppC_v2csreprss_reduced_interaction.tex")


# UN reprisals
f0 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1 * v2x_polyarchy_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1 * v2x_polyarchy_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1 * v2x_polyarchy_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1 * v2x_polyarchy_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')


##### Table A.16 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3
                          , type="latex"
                          , title = "Reprisal against UN collaborators and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)) ,
                                             c("Country fixed effects", "No", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "No", "No", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes"))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Polyarchy (1 yr lag)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP communications * Polyarchy"
                                              , "Constant")                          , column.labels = c("(1)","(2)","(3)","(4)", "(5)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Reprisal (original data)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:reprisals_reduced_interaction"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{13cm}{ \\textit{Notes:} The models in columns 1, 2 and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_reduced_interaction.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA16_AppC_unreprisal_acc_reduced_interaction.tex")


### Additional model specifications with different control variables
controls <- paste0("+ v2peedueq_lag1 + conflict_lag + v2x_polyarchy + log_popsize_lag1 + theta_mean_lag1")

# Models of repression of CSOs

f0 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 | 0 | 0 | cowc"))
f1 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("v2csreprss ~  commcount_individual_lag1", controls, "+ lag(v2csreprss,1) | lag(v2csreprss,2:99)") )
f4 <- as.formula( paste0("v2csreprss ~  1 ", controls , " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1) | cowc"))

m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data)

# First stage (save for analysis_V3_first_stage.R to make the table with all first stage regression models)
first <- felm(as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, " + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4) )
summary(first)

# Summary statistics
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.3 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Repression of CSOs and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Polyarchy (lag 1 yr)"
                                              , "Logged population size (lag 1 yr)"
                                              , "Fariss HR scores"
                                              , "Repression of CSOs (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications (lag 1 yr,  instr.)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Repression of CSOs (V-Dem)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:v2csreprss_full"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. \\ $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/v2csreprss_full_altmodels.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA3_AppC_v2csreprss_full_altmodels.tex")


# Models of reprisals against UN collaborators
f0 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1 | cowc | 0 | cowc"))
f1 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1 | cowc | 0 | cowc"))
f2 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1", controls,  " | cowc | 0 | cowc"))
f3 <- as.formula( paste0("unreprisal_acc ~ commcount_individual_lag1", controls, "+ lag(unreprisal_acc,1) | lag(unreprisal_acc,2:99)") )
f4 <- as.formula( paste0("unreprisal_acc ~ ", controls, " + spvisits_lag1 | cowc | (commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1 ) | cowc"))

# Model estimation
m0 <- felm(formula = f0, data = data )
m1 <- felm(formula = f1, data = data )
m2 <- felm(formula = f2, data = data )
m3 <- pgmm(formula = f3, data = data, index=c("cowc", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = data )

# Summary stats (note that first stage is same as above)
first <- felm( as.formula( paste0("commcount_individual_lag1 ~ new_entries_sum_lag1 + Rapporteurs_active_number_lag1", controls, "  + spvisits_lag1 | cowc | 0 | cowc")), data = model.frame(m4))
summary(first)
fstat <- round(summary(first)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  

##### Table A.4 for Appendix C #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , title = "Reprisal for cooperation with UN and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", " ", " ", " ", " ",  fstat))
                       , covariate.labels = c("UNSP communications (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Polyarchy (lag 1 yr)"
                                              , "Logged population size (lag 1 yr)"
                                              , "Fariss HR score"
                                              , "Reprisal for cooperation with UN (lag 1 yr, instr.)"
                                              , "UNSP country visits (lag 1 yr)"
                                              , "UNSP communications (lag 1 yr,  instr.)")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Reprisal against UN collaborators (original data)"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:unreprisal_acc_full"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/unreprisal_acc_full_altmodels.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA4_AppC_unreprisal_acc_full_altmodels.tex")


### In which country-years were there no communications?
dat_comm <- data %>% select(cowc, year, cowc_year, commcount_individual) %>% filter(!is.na(commcount_individual)) 
dat_comm_zero <- data %>% select(cowc, year, cowc_year, commcount_individual) %>% filter(!is.na(commcount_individual)) %>% filter(commcount_individual == 0)
nrow(dat_comm_zero) /  nrow(dat_comm)
### In which country there was no communication over the entire analysis period?
least_comm <- dat_comm_zero %>% count(cowc) %>% arrange(desc(n))



#### Check endogeneity by exchange IV and DV #####

f_endo1 <- as.formula( paste0("commcount_individual ~  v2csreprss_lag1", controls,  " | cowc | 0 | cowc"))
m_endo1 <- felm(formula = f_endo1 , data = data )
summary(m_endo1)

f_comp1 <- as.formula( paste0("v2csreprss ~  unreprisal_acc_lag1", controls,  " | cowc | 0 | cowc"))
m_comp1 <- felm(formula = f_comp1 , data = data )
summary(m_comp1)

f_endo2 <- as.formula( paste0("commcount_individual ~  unreprisal_acc_lag1", controls,  " | cowc | 0 | cowc"))
m_endo2 <- felm(formula = f_endo2 , data = data )
summary(m_endo2)

f_comp2 <- as.formula( paste0("unreprisal_acc ~  commcount_individual_lag1", controls,  " | cowc | 0 | cowc"))
m_comp2 <- felm(formula = f_comp2 , data = data )
summary(m_comp2)


##### Table A.17 for Appendix C #####
star_full <- stargazer(m_endo1, m_endo2
                       , title = "Models of UNSP communications as function of CSO repression and reprisal"
                       , digits=3
                       , font.size="scriptsize"
                       , omit.stat = c("n", "ser")
                       , add.lines = list(c("Number of observations"
                                            , length(m_endo1$fitted.values)
                                            , length(m_endo1$fitted.values) ) ,
                                          c("Country fixed effects", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "Yes"),
                                          c("Lagged depend. var.", "No", "No"))
                       , covariate.labels = c("Repression against CSOs (lag 1 yr)"
                                              , "Reprisals against UN collaborators (lag 1 yr)"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "Logged GDP per capita (lag 1 yr)")
                       , column.labels = c("(1)","(2)" )
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Number of UNSP communications"
                       , column.sep.width = "1pt"
                       , align=F
                       , label = "tab:endogeneity"
                       , model.numbers=F, model.names=F
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{3}{l} {\\parbox[t]{15cm}{ \\textit{Notes:} The models in columns 1 and 2 are estimated with ordinary least squares regression. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/endogeneity_models_of_communications.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA17_AppC_endogeneity_models_of_communications.tex")



### exit ####


